!-------------------------------------------------------------------------------------
subroutine diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, &          ! inputs 
                        t, Courant_in_z, OmegaBounce)            ! inputs
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : kz, pi, nz, runname, restart, SpeckleWidth
  use efield, only     : efield_get_epsilon, efield_get_E0z_center

  implicit none

  !-- global variables --

  real, intent(in)       :: dfu(jv, jz, jy, jx, jDir)
  real, intent(in)       :: Ez2D(jz,jy,jx), Ex2D(jz,jy,jx), Ey2D(jz,jy,jx)
  real, intent(in)       :: f0vec(jv)
  real, intent(in)       :: t, Courant_in_z, OmegaBounce


  !-- local variables --

  logical, save          :: first_call = .true.

  real       :: PhasePhi
  real, save :: PhasePhiOld
  real       :: KE, PotE, PotEx, PotEy, potEz, Phi0Energy
  real       :: nTrapped, fsq, fsq1


  !complex    :: epsilon
  complex    :: EzFtInternal, EzFTtotal

            
  real       :: funiform(jv), funiform1(jv)

  complex    :: Phi1Dz(jz)
  real       :: E0zCenter(jz)

  integer    :: ix, iy, iz, iv, iDir  ! counters


  !-----------------------------------------------------------------------------------
       
       
      fsq1 = 0  
      do iv=1,jv
        funiform1(iv) = sum( dfu(iv,:,:,:,:) )
        fsq1 =  fsq1 + sum( ( f0vec(iv) + dfu(iv,:,:,:,:) )**2  )
      end do

      call mpi_reduce(funiform1, funiform, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      call mpi_reduce(fsq1, fsq, 1, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)


      if (myid .ne. 0) return
     
      funiform = funiform /(jx*jy*jz)

      fsq = fsq*dz*dv 
      if (jx .ne. 1) fsq = fsq*dx
      if (jy .ne. 1) fsq = fsq*dy


      PotEx = 0.5*sum(Ex2D(:,:,:)**2) / (jx*jy*jz)
      PotEy = 0.5*sum(Ey2D(:,:,:)**2) / (jx*jy*jz)       
      PotEz = 0.5*sum(Ez2D(:,:,:)**2) / (jx*jy*jz)

      potE = PotEx + PotEy + PotEz
        

      KE =  0.5 * sum( funiform(:) * vvec(:)**2 ) * dv

      call efield_get_E0z_center(E0zCenter)

      Phi1Dz(:) = Ez2D(:, jy/2+1, jx/2+1) - E0zCenter(:)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFtInternal = Phi1Dz(1+nz)
        
      Phi1Dz(:) = Ez2D(:, jy/2+1, jx/2+1)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFTtotal = Phi1Dz(1+nz)

      PhasePhi = -aimag(log(EzFtInternal)) ! really phase of E
        
      do while ( ( PhasePhi - PhasePhiOld ) <  - pi )
         PhasePhi = PhasePhi + 2*pi
      enddo
      do while ( ( PhasePhi - PhasePhiOld ) >  + pi )
         PhasePhi = PhasePhi - 2*pi
      enddo
        
      PhasePhiOld = PhasePhi   

     !------------------------------
       

      if (first_call .and. .not. restart) then

         open(3, file=trim(runname)//'.dat')

         write(3, "('% 1.t  2.PhiInternal  3.PhiToTal  4.-Phase(phi)  5.KineticEnergy')")
         write(3, "('% 6.Ex^2energy  7.Ey^2energy  8.Ez^2energy  9.OmegaBounce  10.Courant_in_z')")
         write(3,*)

      else

         open(3, file=trim(runname)//'.dat', position='append')

      end if

      if (first_call)  first_call = .false.

      PHi0Energy=0
      nTrapped=0.0

           write(3,501) t, 2*abs(EzFtInternal)/kz, 2*abs(EzFTtotal)/kz, &
           PhasePhi, KE, potEx, potEy, potEz, &
           OmegaBounce*SpeckleWidth*0.37233, Courant_in_z

      close(3)

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics


!--------------------------------------------------------------------------------------
subroutine write_modulation(xMod, uxMod, yMod, uyMod)
!--------------------------------------------------------------------------------------

  use dimensions
  use parameters, only :  runname

  implicit none

  complex, intent(in)  :: xMod(jx), uxMod(jx), yMod(jy), uyMod(jy)
  integer              :: ix, iy

  open(1, file = trim(runname)//'.modulation.txt')

  write(1,*) '% 1.x    2.xMod   3.uxMod'
  write(1,*)

  do ix =1, jx
     write(1,501) xvec(ix), real(xMod(ix)),  real(uxMod(ix))
  enddo

  write(1,*)
  write(1,*)


  write(1,*) '% 1.x    2.yMod   3.uyMod'
  write(1,*)

  do iy =1, jy
     write(1,501) yvec(iy), real(yMod(iy)),  real(uyMod(iy))
  enddo

  close(1)


501   format(  6400(e12.6,2x)  )


end subroutine write_modulation

!--------------------------------------------------------------------------------------
subroutine write_dispersion()
!--------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap
  use parameters, only :  kz, runname

  implicit none

  integer              :: iz
  real                 :: phi0, phi00
  real                 :: Omega0, DeltaOmega, dOmegadPhi


  Omega0(kz) = 1.106 + 1.969*(kz-0.25) +3.16*(kz-0.25)**2 - 13.0*(kz-0.25)**3  & 
       + 12.11*(kz-0.25)**4 + 1.068*(   (0.534-kz)**0.5  -  (0.534-0.25)**0.5  )
  ! accurate for 0.25<kz<0.53
 
  DeltaOmega(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * &
       ( exp(-0.5 / kz**2) / kz**2 ) * phi0**0.5  & 
       - 343 * kz**4 * exp( -0.5 / kz**2 ) * phi0**1.5
  !  for   0.275 < kz < 0.425
 
  dOmegadPhi(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * &
       ( exp(-0.5 / kz**2) / kz**2 ) * 0.5 / phi0**0.5  & 
       - 343 * kz**4 * exp( -0.5 / kz**2 ) * 1.5 * phi0**0.5
  !  for   0.275 < kz < 0.425


  if (myid .ne. 0) return
    
  open(1, file = trim(runname)//'.dispersion.txt')
  write(1,*)'% 1.phi   2.Vphase     3.-0.5*phixdomega/dphi'
  write(1,*)
  do iz=1,40
     phi00 = iz * 0.5/40
     write(1,501) phi00, ( Omega0(kz) + DeltaOmega(kz,phi00 ) ) /  kz, & 
          -0.5 * phi00 * dOmegadPhi(kz,phi00)
  end do
  close (1)


501   format(  6400(e12.6,2x)  )

end subroutine write_dispersion

!-------------------------------------------------------------------------------------
subroutine write_f0_initial(f0vec, dfu)
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap
  use parameters, only   : runname

  implicit none

  real, intent(in)       :: dfu(jv, jz, jy, jx, jDir)
  real, intent(in)       :: f0vec(jv)

  real                   :: meandf1(jv), meandf(jv)
  integer                :: iv
 

  do iv=1,jv
     meandf1(iv) = sum(dfu(iv,:,:,:,:)) / (jx*jy*jz)
  enddo

  call mpi_reduce(meandf1, meandf, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

  if (myid .ne. 0) return
 
  open(13, file = trim(runname)//'.f0initial.txt')
  write(13,*) '% 1.V  2.f0  3.<df>'
  write(13,*) 

  do iv=1,jv
     write(13,501) vvec(iv), f0vec(iv),  meandf(iv)
  enddo

  close(13)

501   format(  6400(e12.6,2x)  )

end subroutine write_f0_initial

!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------


